Df <- fread("N:/durable/projects/crayner/PipePsy/Data/ChPsyFsc.tsv") %>% select(PrID,starts_with("Ch"))
names(Df) <- str_replace_all(names(Df),"Ch","")
names(Df) <- str_replace_all(names(Df),"Q","_Q")
names(Df) <- str_replace_all(names(Df),"Q04","6mo")
names(Df) <- str_replace_all(names(Df),"Q05","18mo")
names(Df) <- str_replace_all(names(Df),"Q06","3yr")
names(Df) <- str_replace_all(names(Df),"Q07","5yr")
names(Df) <- str_replace_all(names(Df),"Q09","8yr")
names(Df) <- str_replace_all(names(Df),"_Q11C","C_14yr")
names(Df) <- str_replace_all(names(Df),"_Q11M","M_14yr")
Df <- Df %>% select(PrID,matches("6mo"),matches("18mo"),matches("3yr"),matches("5yr"),matches("8yr"),matches("14yr"))
excl <- c(
"IcqFull_6mo","CbclFull_18mo","EasFull_18mo","CbclFull_3yr","CfqFull_3yr",
"EasFull_3yr","ItseaFull_3yr","ScqFull_3yr","SdqFull_3yr","CbclFull_5yr",
"CccFull_5yr","EasFull_5yr","SprakFull_5yr","CccFull_8yr","NhpicFull_8yr",
"RsdbdFull_8yr","ScqFull_8yr","SprakFull_8yr","PeqFullC_14yr","SbFullC_14yr",
"SclFullC_14yr","SdqFullC_14yr","SdqFullM_14yr","SppaFullC_14yr",
"CbclAdhd_18mo","CbclAdhd_3yr","CbclAdhd_5yr"
)
Df <- Df %>% select(!c(all_of(excl)))
Cov <-
# data.table::fread("N:/durable/projects/crayner/Data/Participation/MoBaAgeVarsLong.csv")%>%
data.table::fread("N:/durable/projects/crayner/Data/Participation/MoBaAgeVarsLongLinked.csv") %>%
select(MoLNR,FaLNR,PrID,ChIID,ChN,ChAge,MoAge,FaAge,Wave=Q) %>% # ChSex=ChSEX,
filter(ChAge>0&Wave!="QF2"&Wave!=""&Wave!=" ") %>%
distinct() %>% # filter(ChAge>0) %>%
mutate(PrID=paste0(PrID,"_",ChN), Wave=droplevels(as.factor(Wave))) %>%
filter(PrID %in% Df$PrID) %>%
mutate(FamID=as.numeric(as.factor(paste0(MoLNR,"_",FaLNR)))) %>%
select(PrID,FamID,ChIID,Wave,matches("Age")) %>% #,ChSex
mutate(Wave=factor(Wave,levels=c("Q04","Q05","Q06","Q07","Q09","Q11C","Q11M"),
labels=c("6mo","18mo","3yr","5yr","8yr","C14yr","M14yr"))) %>%
filter(ChIID!=""&ChIID!=" "&!is.na(ChIID))
AgeDiff <- Cov %>% filter(Wave=="6mo") %>% mutate(ChMo=MoAge-ChAge,ChFa=FaAge-ChAge) %>% select(PrID,ChMo,ChFa)
Cov2 <- data.table::fread("N:/durable/projects/crayner/Data/Covariates/ChAge14yrLong.tsv") %>%
inner_join(AgeDiff) %>%
mutate(MoAge=ChMo+ChAge,FaAge=ChFa+ChAge)%>%select(-ChMo,-ChFa)
Cov3 <- data.table::fread("N:/durable/projects/crayner/Data/Covariates/MbrnCovariates113536.tsv")%>%
mutate(PrID=paste0(PrID,"_",BaN)) %>%
filter(PrID %in% Df$PrID) %>%
mutate(ChSex=factor(ChildSGender,levels=c("Male","Female"),labels=c(1,2))) %>%
select(PrID,ChSex)
Cov <- Cov %>% bind_rows(Cov2) %>% inner_join(Cov3)
tmpDf_1 <-
Df %>%
select(PrID,matches("6mo")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="6mo")) %>%
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")
tmpDf_2 <-
Df %>%
select(PrID,matches("18mo")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.35) %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="18mo")) %>%
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")
tmpDf_3 <-
Df %>%
select(PrID,matches("3yr")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="3yr")) %>%
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")
tmpDf_4 <-
Df %>%
select(PrID,matches("5yr")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="5yr")) %>%
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")
tmpDf_5 <-
Df %>%
select(PrID,matches("8yr")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.1) %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="8yr")) %>%
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t()%>% as.data.frame()%>% rownames_to_column("labs")
tmpDf_6 <-
Df %>%
select(PrID,matches("C_14yr")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% filter(nNA>0.5) %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="C14yr")) %>% #) %>% distinct(PrID,.keep_all=T)
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t() %>% as.data.frame() %>% rownames_to_column("labs")
tmpDf_7 <-
Df %>%
select(PrID,matches("M_14yr")) %>%
mutate(nNA=apply(X=!is.na(.),MARGIN=1,FUN=mean)) %>% na.omit() %>% #na.omit() %>%
select(PrID) %>% inner_join(Cov %>% filter(Wave=="M14yr")) %>% #) %>% distinct(PrID,.keep_all=T)
summarize( N=n(),
`Males`=sum(ChSex==1,na.rm=T),
`Males %`=sum(ChSex==1,na.rm=T)/n(),
`Females`=sum(ChSex==2,na.rm=T),
`Females %`=sum(ChSex==2,na.rm=T)/n(),
`Child Age (mean)`=round(mean(ChAge),2),
`Mothers Age (mean)`=round(mean(MoAge),2),
`Fathers Age (mean)`=round(mean(FaAge),2),
`Child Age (sd)`=round(sd(ChAge),2),
`Mothers Age (sd)`=round(sd(MoAge),2),
`Fathers Age (sd)`=round(sd(FaAge),2)) %>%
select(N,`Males`,`Females`,`Males %`,`Females %`,matches("Ch"),matches("Mo"),matches("Fa")) %>% t() %>% as.data.frame() %>% rownames_to_column("labs")
tmpDf <- plyr::join_all(list(tmpDf_1,tmpDf_2,tmpDf_3,tmpDf_4,tmpDf_5,tmpDf_6,tmpDf_7),by="labs") %>% column_to_rownames("labs")
names(tmpDf)<- c("6mo","18mo","3yr","5yr","8yr","C14yr","M14yr")
datatable(
tmpDf %>% mutate(across(everything(), function(x) round(x,2))),
extensions='Buttons',
options=list(
dom='tiB',buttons=c('csv'),
searching=F,paging=F,ordering=F,
columnDefs=list(list(className='dt-left',targets="_all"))
))
Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.
Df <- fread("N:/durable/projects/crayner/PipePsy/Data/ChPsyFsc.tsv") %>% select(PrID,starts_with("Ch"))
names(Df) <- str_replace_all(names(Df),"Ch","")
names(Df) <- str_replace_all(names(Df),"Q","_Q")
names(Df) <- str_replace_all(names(Df),"Q04","6mo")
names(Df) <- str_replace_all(names(Df),"Q05","18mo")
names(Df) <- str_replace_all(names(Df),"Q06","3yr")
names(Df) <- str_replace_all(names(Df),"Q07","5yr")
names(Df) <- str_replace_all(names(Df),"Q09","8yr")
names(Df) <- str_replace_all(names(Df),"_Q11C","C_14yr")
names(Df) <- str_replace_all(names(Df),"_Q11M","M_14yr")
Df <- Df %>% select(PrID,matches("6mo"),matches("18mo"),matches("3yr"),matches("5yr"),matches("8yr"),matches("14yr"))
excl <- c(
"IcqFull_6mo","CbclFull_18mo","EasFull_18mo","CbclFull_3yr","CfqFull_3yr",
"EasFull_3yr","ItseaFull_3yr","ScqFull_3yr","SdqFull_3yr","CbclFull_5yr",
"CccFull_5yr","EasFull_5yr","SprakFull_5yr","CccFull_8yr","NhpicFull_8yr",
"RsdbdFull_8yr","ScqFull_8yr","SprakFull_8yr","PeqFullC_14yr","SbFullC_14yr",
"SclFullC_14yr","SdqFullC_14yr","SdqFullM_14yr","SppaFullC_14yr",
"CbclAdhd_18mo","CbclAdhd_3yr","CbclAdhd_5yr"
)
Df <- Df %>% select(!c(all_of(excl)))
# excl <- names(Df)[grepl(x=names(Df),"Full")]
# fullScales <- c(
# "CastFull_5yr","CprsrsFull_5yr","SlasFull_5yr","CebqFull_8yr","ScaredFull_8yr",
# "SmfqFull_8yr","SpinFull_8yr","CapeFullC_14yr","SmfqFullC_14yr","SwlsFullC_14yr",
# "YptiFullC_14yr"
# )
# fullScales <- c("Cast","Cprsrs","Slas","Cebq","Scared","Smfq","Spin","Cape","Smfq","Swls","Ypti")
# test <- names(Df)[grepl(x=names(Df),paste0(fullScales,collapse="|"))]
Tab1 <- lapply(Df, function(x) sum(!is.na(x)))
Tab1 <- do.call(Tab1,what="bind_rows")
Tab1 <- t(Tab1) %>% as.data.frame()
Tab1 <- Tab1 %>% rownames_to_column("Var")
Tab1 <- Tab1[-1,]
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"C_14yr", "_14yr Child"))
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"M_14yr", "_14yr Mother"))
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"([a-z])([A-Z])", "\\1 \\2"))
Tab1 <- Tab1 %>% mutate(Var=str_replace_all(Var,"_", " "))
Tab1 <- Tab1 %>% separate(Var, c("Scale","Subscale","Age","Rater"))
Tab1 <- Tab1 %>% mutate(Rater=ifelse(is.na(Rater),"Mother",Rater))
Tab1 <- Tab1 %>% mutate(Scale=str_to_upper(Scale))
datatable(
Tab1 %>% select(Scale,Subscale,Age,Rater,N=V1),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
ICQ6, Infant Characteristics Questionnaire 6 Months Form; ASQ, Ages and Stages Questionnaires; EAS, Emotionality, Activity and Shyness Temperament Questionnaire; CBCL, Child Behaviour CheckList; AUTQ, Early Screening of Autistic Traits Questionnaire-Short (+Modified Checklist for Autism in Toddlers); ITSEA, Infant-Toddler Social and Emotional Assessment; CFQ, The Child Feeding Questionnaire; SDQ, Strength and Difficulties Questionnaire-Prosocial Subscale; SCQ, Social Communication Questionnaire; PPBS, Preschool Play Behaviour Scale; CPRSRS, Conners Parent Rating Scale - Revised, Short Form; CCC, Children’s Communication Checklist-2 Short: Coherence I; SPRAK, Statements about Language-Related Difficulties: Semantics; SLAS, Speech and Language Assessment Scale; CAST, Childhood Asperger Syndrome Test; RSDBD, Rating Scale for Disruptive Behaviour Disorders; CEBQ, Children’s Eating Behaviour Questionnaire; SCARED, Screen for Child Anxiety Related Disorders; SMFQ, Short Mood and Feelings Questionnaire; SPIN, Mini Social Phobia Inventory; NHPIC30, Short Norwegian Hierarchical Personality Inventory for Children; PEQ, Parental Environment Questionnaire; YPTI, Youth Psychopathic Traits; CAPE15, Community Assessment of Psychic Experiences; SCL, The (Hopkins) Symptoms Checklist; DES, Differential Emotional Scale; RSES, Rosenberg Self-Esteem Scale; SWLS, Satisfaction with life scale; IPIP, International Personality Item Pool; SFS, School Functioning Scale; SPPA, Self-Perception Profile for Adolescents; SB, School belonging
ModStat <- fread("Tables_Phe/xIRT_ModelStats.tsv") %>%
filter(grepl("Ch",Scale)) %>%
mutate(Trait=paste0(Scale,Subscale)) %>%
select(Trait,Age=Wave,everything()) %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
mutate(Age=str_replace_all(Age,"Q04","6mo")) %>%
mutate(Age=str_replace_all(Age,"Q05","18mo")) %>%
mutate(Age=str_replace_all(Age,"Q06","3yr")) %>%
mutate(Age=str_replace_all(Age,"Q07","5yr")) %>%
mutate(Age=str_replace_all(Age,"Q09","8yr")) %>%
mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>%
mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>%
select(!matches("Scale"),!matches("Subscale")) %>%
filter(!Trait %in% excl) %>%
select(-Scale,-p,-RMSEA_5,-RMSEA_95,-LL)
datatable(
ModStat,
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=T,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,scrollY=T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
1 See table one for description of each trait. 2 Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c. 3 The item response theory model which had the best fit to the data: graded, graded response model (ref); nominal, nominal response model; 3-PL, three parameter logistic model (ref); Rasch model (ref). 4AIC, Akaike Information Criterion. 5 BIC, Bayesian Information Criterion. 6 Model fit assessed using M2 index (2 variant), which is specifically designed to assess the fit of item response models for ordinal data - when polytomous response models do not have sufficient degrees of freedom to compute M2. This function also computes associated fit indices that are based on fitting the null model. 7 The degrees of freedom for the M2 test. 8 RMSEA, the M2-based root mean square error of approximation is the primary fit index - can be used to suggest that the data fit the model, with suggested cutoff values of RMSEA <= .06. 9 SRMSR, standardised root mean square residual and 10 CFI, comparative fit index is used to assess adequacy of model fit. The suggested cutoff values of SRMSR <= .08. 12 Reliability, a measure of the average inter-item correlation, assessed using Cronbach’s alpha.
# ESTIMATING MATRICES (HETCOR)
CorMat <- cormat(Df %>% select(-PrID))
# Make interactive correlation table
CorTably <- cortab(CorMat)
# Plot interactive correlation matrix
CorTab <- cortably(CorMat)
CorTab$z <- abs(CorTab$r/CorTab$se)
Cor <-
CorTab %>%
select(Var1,Var2,r) %>%
filter(Var1!=Var2) %>%
pivot_wider(names_from = Var2, values_from = r) %>%
column_to_rownames("Var1")
Err <-
CorTab %>%
select(Var1,Var2,se) %>%
filter(Var1!=Var2) %>%
pivot_wider(names_from = Var2, values_from = se) %>%
column_to_rownames("Var1")
Zsc <-
CorTab %>%
select(Var1,Var2,z) %>%
filter(Var1!=Var2) %>%
pivot_wider(names_from = Var2, values_from = z) %>%
column_to_rownames("Var1")
Zsc <- log10(Zsc)
ChCorplotlyFull <- heatmaply::heatmaply_cor(
as.data.frame(Cor),
node_type = "scatter",
point_size_mat = Zsc,
point_size_name = "log10|Z|",
label_names = c("x", "y", "r"),
showticklabels = c(T,F),
fontsize_col=7,
grid_size = 0.001,
dendrogram = "both",
subplot_heights=c(0.04, 0.96),
subplot_widths=c(0.94,0.06),
hide_colorbar = T
)
ChCorplotlyFull$height <- nrow(Cor)*13
ChCorplotlyFull$width <- nrow(Cor)*13
# ChCorplotlyFull$x$data[[2]]$marker$size <- (ChCorplotlyFull$x$data[[2]]$marker$size/3)
# ChCorplotlyFull$x$data[[3]]$marker$size <- (ChCorplotlyFull$x$data[[3]]$marker$size/3)
ChCorplotlyFull$x$data[[4]]$marker$size <- (ChCorplotlyFull$x$data[[4]]$marker$size/2.5)
# ChCorplotlyFull$x$data[[5]]$marker$size <- (ChCorplotlyFull$x$data[[5]]$marker$size/3)
ChCorplotlyFull
vDf <- fread("N:/durable/projects/crayner/PipePsy/Data/FamPsyQrs.tsv") %>% select(PrID,starts_with("Ch"))
names(vDf) <- str_replace_all(names(vDf),"Ch","")
names(vDf) <- str_replace_all(names(vDf),"Q","_Q")
names(vDf) <- str_replace_all(names(vDf),"Q04","6mo")
names(vDf) <- str_replace_all(names(vDf),"Q05","18mo")
names(vDf) <- str_replace_all(names(vDf),"Q06","3yr")
names(vDf) <- str_replace_all(names(vDf),"Q07","5yr")
names(vDf) <- str_replace_all(names(vDf),"Q09","8yr")
names(vDf) <- str_replace_all(names(vDf),"_Q11C","C_14yr")
names(vDf) <- str_replace_all(names(vDf),"_Q11M","M_14yr")
vDf <- vDf %>% select(PrID,matches("6mo"),matches("18mo"),matches("3yr"),matches("5yr"),matches("8yr"),matches("14yr"))
vDf <- vDf %>% select(!c(excl))
# ESTIMATING MATRICES (HETCOR)
# vCorMat <- cormat(vCorvDf)
vCorMat <- cormat(vDf %>% select(-PrID))
# Make interactive correlation table
vCorTably <- cortab(vCorMat)
# Plot interactive correlation matrix
vCorTab <- cortably(vCorMat)
vCorTab$z <- abs(vCorTab$r/vCorTab$se)
vCor <-
vCorTab %>%
select(Var1,Var2,r) %>%
filter(Var1!=Var2) %>%
pivot_wider(names_from = Var2, values_from = r) %>%
column_to_rownames("Var1")
vCor[is.na(vCor)] <- 0
Err <-
vCorTab %>%
select(Var1,Var2,se) %>%
filter(Var1!=Var2) %>%
pivot_wider(names_from = Var2, values_from = se) %>%
column_to_rownames("Var1")
Err[is.na(Err)] <- 0
Zsc <-
vCorTab %>%
select(Var1,Var2,z) %>%
filter(Var1!=Var2) %>%
pivot_wider(names_from = Var2, values_from = z) %>%
column_to_rownames("Var1")
Zsc <- log10(Zsc)
#TO DO: ADD SE LABELS
ChvCorplotlyFull <- heatmaply::heatmaply_cor(
as.data.frame(vCor),
node_type = "scatter",
point_size_mat = Zsc,
point_size_name = "log10|Z|",
# point_size_mat = logP,
# point_size_name = "-log10P/100",
label_names = c("x", "y", "r"),
showticklabels = c(T,F),
fontsize_col=7,
grid_size = 0.001,
dendrogram = "both",
subplot_heights=c(0.06, 0.94),
subplot_widths=c(0.94,0.06),
hide_colorbar = T
)
ChvCorplotlyFull$height <- nrow(vCor)*13
ChvCorplotlyFull$width <- nrow(vCor)*13
# ChvCorplotlyFull$x$data[[2]]$marker$size <- (ChvCorplotlyFull$x$data[[2]]$marker$size/3)
# ChvCorplotlyFull$x$data[[3]]$marker$size <- (ChvCorplotlyFull$x$data[[3]]$marker$size/3)
ChvCorplotlyFull$x$data[[4]]$marker$size <- (ChvCorplotlyFull$x$data[[4]]$marker$size/2.5)
# ChvCorplotlyFull$x$data[[5]]$marker$size <- (ChvCorplotlyFull$x$data[[5]]$marker$size/3)
ChvCorplotlyFull
screePlot <- scree_plot(data=Df%>%select(-PrID), method = "ml", n.iter=200)+ scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot1 <- screePlot + scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot1
screePlot <- scree_plot(data=vDf%>%select(-PrID), method = "ml", n.iter=200)+ scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot2 <- screePlot + scale_x_continuous(breaks=1:25,limits=c(1,25))
screePlot2
PLOT1
PLOT2
PLOT3
PLOT4
datatable(
CorTably %>% select(Var1,Var2,r,se,p,n),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=T,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,scrollY=T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
datatable(
vCorTably %>% select(Var1,Var2,r,se,p),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
Violin plots display the distribution of a phenotype by minor allele count, at genotypes with different effects on the mean and variance of the phenotype. (A-D) The red dotted lines represent the regression of phenotype on genotype. Blue lines represent regression of phenotype on genotype, restricted to individuals in the upper and lower quintiles of the phenotype distribution - this helps to visualise the heteroscedasticity across genotype levels. All genotypes were simulated with minor allele frequency equal to 0.4 and environments had a normal distribution (n=1000). All phenotypes were simulated as a function of genetic and/or environmental effects with random normal error distribution. (A) For the null model the phenotype was simulated without any effect from genotype or environment. (B) the mQTL model was simulated with genetic effect 𝛽=.2. (C) the vQTL model was simulated with genetic effect 𝛽=0, an environmental effect 𝛽=.5 and a GxE interaction effect 𝛽=.4. (D) the mvQTL model was simulated with genetic effect 𝛽=.2, an environmental effect 𝛽=.5 and a GxE interaction effect 𝛽=.4. (E-F) Points are coloured based on a second variable. Lines represent regression of phenotype on genotype stratified by (E) tertiles of environmental exposure and (F) a genotype at a second locus. This illustrates how the interaction between a genotype and a second factor (an environmental exposure or a second genetic locus) contributes to increased variance across the genotype levels. (E) the GxE model was simulated with genetic effect 𝛽=.2, an environmental effect 𝛽=.5 and a GxE interaction effect 𝛽=.4. (F) the GxG model was simulated with genotype 1 effect 𝛽=.2, genotype 2 effect 𝛽=.2 and a GxG interaction effect 𝛽=.2.
Quantile-quantile plots (left) display the expected and observed p-values from GWA (+) and vGWA (x). Scatter plots (right) show the relationship between p-values for SNPs estimated from vGWA on the y-axis and GWA on the x-axis. The relationship between p-values for ASQ motor skills is strongly positively correlated (i.e. the same SNPs are identified as having effects), whereas for ITSEA emotional dysregulation SNPs identified as having stronger variance effects did not necessarily have mean effects.
SnpTab <- fread("Tables_Snps/ChPsy/ClumpedTopSnps.tsv")
SnpTab$af <- round(SnpTab$af,2)
SnpTab$mP <- formatC(SnpTab$mP, format = "e", digits = 2)
SnpTab$vP <- formatC(SnpTab$vP, format = "e", digits = 2)
DT::datatable(
SnpTab,
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
CorP <- cor.test(as.numeric(SnpTab$mP),as.numeric(SnpTab$vP))
1 See table one for description of each trait. 2 Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.
Correlation between p-values for mGWA and vGWA: Pearson’s r = -0.36; p = 7.16e-01
vPheWasTab <- fread("Tables_Snps/ChPsy/vPheWasTab.tsv")
mPheWasTab <- fread("Tables_Snps/ChPsy/mPheWasTab.tsv")
PheWasTab3 <- vPheWasTab %>% bind_rows(mPheWasTab)
PheWasLab <-
PheWasTab3 %>%
arrange(PheWasP) %>% distinct(PheWasTrait,MoBaTrait,.keep_all=T) %>%
select(PheWasTrait,MoBaTrait,PheWasP) %>%
mutate(Label=str_replace_all(PheWasTrait, "\\s*\\([^\\)]+\\)", "")) %>%
mutate(Label=str_replace_all(Label, "[[:punct:]]", "")) %>%
mutate(Label=str_to_title(Label)) %>%
mutate(Label=str_replace_all(Label, " ", "")) %>%
mutate(Label=str_replace_all(Label, "CellCount", "Count")) %>%
distinct(Label,MoBaTrait,.keep_all=T) %>%
inner_join(PheWasTab3,c("MoBaTrait","PheWasTrait","PheWasP")) %>%
select(Label,PheWasTrait,MoBaTrait,vP,mP,PheWasP,rsid) %>% distinct() %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "C_", " ")) %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "Full", "")) %>%
# mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([A-Z])", "\\1 \\2")) %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([0-9])", "\\1 \\2"))%>%
separate(MoBaTrait,c("MoBaTrait","Age")," ")%>%
mutate(Age = factor(Age, levels=c("6mo","18mo","3yr","5yr","8yr","14yr"))) #%>%
# mutate(MoBaTrait=paste0(Scale,Trait))
PheWasLab <-
bind_rows(
PheWasLab %>% filter(vP<5e-8) %>% rename("GwaP"=vP) %>% mutate(Effect="vQTL"),
PheWasLab %>% filter(mP<5e-8) %>% rename("GwaP"=mP) %>% mutate(Effect="mQTL")
) %>%
mutate(GwaP=round(-log10(GwaP),1),PheWasP=round(-log10(PheWasP),1))
PheWasTab <- PheWasTab3 %>% select(rsid,PheWasTrait,MoBaTrait,PheWasP,mP,vP)
PheWasTab$PheWasP <- formatC(PheWasTab$PheWasP, format = "e", digits = 2)
PheWasTab$mP <- formatC(PheWasTab$mP, format = "e", digits = 2)
PheWasTab$vP <- formatC(PheWasTab$vP, format = "e", digits = 2)
PheWasPlot <-
ggplot(
PheWasLab,aes(y=PheWasP,x=GwaP,shape=MoBaTrait,colour=MoBaTrait,label=PheWasTrait))+
geom_point()+theme_minimal()+
facet_grid(cols=vars(Age),rows=vars(Effect))+
scale_colour_viridis_d(begin=.2,end=.8)+
scale_shape_manual(values=c(0:14))+
theme(legend.position="bottom")
library(plotly)
ggplotly(PheWasPlot, tooltip = c("MoBaTrait", "GwaP", "PheWasTrait", "PheWasP"))
# PheWasPlotly <- plotly::ggplotly(PheWasPlot)
# saveWidget(PheWasPlotly, "PsyPheWasPlotly.html", selfcontained = T, libdir = "lib")
Scatter plots display the -log10 GWA p-values from publicly available datasets (PheWasP; y-axis) and -log10 GWA p-values from MoBa analyses (x-axis). The plot is stratified by Age of MoBa participants, with age increasing from left to right (x-strip) and by the type of analysis (y-strip), with mGWA in the top panel and vGWA in the bottom panel. MoBa traits are labelled by the shape and the colour of the point. PheWas traits are labelled in hover over boxes (interactive version) and can be viewed in S.Table 5.
DT::datatable(
PheWasTab,
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
1 See table one for description of each trait. 2 Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.
h2TabM <-
fread(paste0("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenFscAsr_h2FullTable.txt"),header=T) %>%
distinct(Trait,.keep_all = T) %>%
mutate(Trait=str_replace_all(Trait,"Q0","_Q0")) %>%
mutate(Trait=str_replace_all(Trait,"Q1","_Q1")) %>%
separate(Trait,c("Trait","Age"),"_") %>%
mutate(alpha = ifelse(h2<0 | h2-SE<0, 1, 1)) %>%
mutate(Type = paste("GWA")) %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
mutate(Age=str_replace_all(Age,"Q04","6mo")) %>%
mutate(Age=str_replace_all(Age,"Q05","18mo")) %>%
mutate(Age=str_replace_all(Age,"Q06","3yr")) %>%
mutate(Age=str_replace_all(Age,"Q07","5yr")) %>%
mutate(Age=str_replace_all(Age,"Q09","8yr")) %>%
mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>%
mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>%
filter(!Trait %in% excl)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
h2TabV <-
fread(paste0("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenQrsAsr_h2FullTable.txt"),header=T) %>%
distinct(Trait,.keep_all = T) %>%
mutate(Trait=str_replace_all(Trait,"Q0","_Q0")) %>%
mutate(Trait=str_replace_all(Trait,"Q1","_Q1")) %>%
separate(Trait,c("Trait","Age")) %>%
mutate(alpha = ifelse(h2<0 | h2-SE<0, 1, 1))%>%
mutate(Type = paste("vGWA")) %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
mutate(Age=str_replace_all(Age,"Q04","6mo")) %>%
mutate(Age=str_replace_all(Age,"Q05","18mo")) %>%
mutate(Age=str_replace_all(Age,"Q06","3yr")) %>%
mutate(Age=str_replace_all(Age,"Q07","5yr")) %>%
mutate(Age=str_replace_all(Age,"Q09","8yr")) %>%
mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>%
mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>%
filter(!Trait %in% excl)
h2Tab1 <- rbind(h2TabM,h2TabV)
h2Tab2 <-
h2TabM %>%
select(Trait,Age,h2_mGWA=h2,SE_mGWA=SE,Z_mGWA=Z, I_mGWA=Int, X2_mGWA=MeanChi2, L_mGWA=Lambda) %>%
full_join(
h2TabV %>%
select(Trait,Age,h2_vGWA=h2,SE_vGWA=SE,Z_vGWA=Z, I_vGWA=Int, X2_vGWA=MeanChi2, L_vGWA=Lambda)
)
N <- fread("Tables_GenSem/ChPsy/ChPsyGenQrsAsr_SampleSizes.txt",header=F) %>%
mutate(V1=str_replace_all(V1,"_RegenieMQt.Munged.sumstats.gz","")) %>%
mutate(V1=str_replace_all(V1,"_RegenieVQt.Munged.sumstats.gz","")) %>%
mutate(V1=str_replace_all(V1,"GwaSumStats/ChPsyGenQrsAsr/","")) %>%
mutate(V1=str_replace_all(V1,"Q", "_Q")) %>%
separate(V1, c("Trait","Age")) %>%
select(Trait, Age, N=V2) %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
mutate(Age=str_replace_all(Age,"Q04","6mo")) %>%
mutate(Age=str_replace_all(Age,"Q05","18mo")) %>%
mutate(Age=str_replace_all(Age,"Q06","3yr")) %>%
mutate(Age=str_replace_all(Age,"Q07","5yr")) %>%
mutate(Age=str_replace_all(Age,"Q09","8yr")) %>%
mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>%
mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>%
filter(!Trait %in% excl)
h2_TabFull <- fread("Tables_Phe/xIRT_ModelStats.tsv") %>%
mutate(Trait=paste0(Scale,Subscale)) %>%
select(Trait,Age=Wave,everything()) %>%
mutate(Rater=ifelse(grepl(x=Age,"C"),"Child","Mother")) %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
mutate(Age=str_replace_all(Age,"Q04","6mo")) %>%
mutate(Age=str_replace_all(Age,"Q05","18mo")) %>%
mutate(Age=str_replace_all(Age,"Q06","3yr")) %>%
mutate(Age=str_replace_all(Age,"Q07","5yr")) %>%
mutate(Age=str_replace_all(Age,"Q09","8yr")) %>%
mutate(Age=str_replace_all(Age,"Q11C","14yr(c)")) %>%
mutate(Age=str_replace_all(Age,"Q11M","14yr(m)")) %>%
select(!matches("Scale"),!matches("Subscale"),!matches("_"),-AIC) %>%
inner_join(h2Tab2) %>% inner_join(N) %>%
filter(!Trait %in% excl) %>%
mutate(Age=factor(Age,levels=c("6mo","18mo","3yr","5yr","8yr","14yr(m)","14yr(c)"),
labels=c("6 m/o","18 m/o","3 y/o","5 y/o","8 y/o","14 y/o","14 y/o")
)) %>%
mutate(
h2_mGWA=round(h2_mGWA,2),h2_vGWA=round(h2_vGWA,2),
SE_mGWA=round(SE_mGWA,2),SE_vGWA=round(SE_vGWA,2),
lSE_vGWA=round(h2_vGWA-SE_vGWA,2),uSE_vGWA=round(h2_vGWA+SE_vGWA,2),
lSE_mGWA=round(h2_mGWA-SE_mGWA,2), uSE_mGWA=round(h2_mGWA+SE_mGWA,2)
)
fwrite(h2_TabFull,"Tables_GenSem/ChPsy/ChPsyh2_FullTable.txt",col.names=T)
meanh2V <- round(mean(h2_TabFull$h2_vGWA),2)
meanZV <- round(mean(h2_TabFull$Z_vGWA),2)
meanh2Vz <- round(mean(h2_TabFull %>% filter(Z_vGWA>3) %>% .$h2_vGWA),2)
meanh2m <- round(mean(h2_TabFull$h2_mGWA),2)
meanZm <- round(mean(h2_TabFull$Z_mGWA),2)
meanh2mz <- round(mean(h2_TabFull %>% filter(Z_mGWA>3) %>% .$h2_mGWA),2)
Figure_3 <-
ggplot(
h2_TabFull,
aes(x=h2_mGWA,y=h2_vGWA)) +
geom_linerange(aes(ymin=round(h2_vGWA-SE_vGWA,2),ymax=round(h2_vGWA+SE_vGWA,2),color=Trait)) +
geom_linerange(aes(xmin=round(h2_mGWA-SE_mGWA,2),xmax=round(h2_mGWA+SE_mGWA,2),color=Trait)) +
scale_x_continuous(limits=c(-0.1,0.2),breaks=c(seq(-.1,.2,0.1)))+
scale_y_continuous(limits=c(-0.1,0.2),breaks=c(seq(-.1,.2,0.1)))+
geom_hline(yintercept=0,colour="limegreen",linetype="dotted") +
geom_vline(xintercept=0,colour="limegreen",linetype="dotted") +
facet_wrap(facets=vars(Age),ncol=3) +
theme_minimal() +
scale_colour_viridis_d(option="plasma") +
guides(color="none")
ggplotly(Figure_3, tooltip = c("Trait","h2_mGWA","SE_mGWA","h2_vGWA","SE_vGWA"))
Scatter plots display the point estimates (intersection) and standard errors for h2mGWA (x-axis) h2vGWA (y-axis) from all MoBa analyses. The plot is stratified by Age of MoBa participants. MoBa traits are labelled in the hover over boxes (interactive version) and can be viewed in S.Table 6.
datatable(
h2_TabFull %>% select(Trait,Age,Rater,h2_mGWA,SE_mGWA,Z_mGWA,h2_vGWA,SE_vGWA,Z_vGWA),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
See table one for description of each trait. 2 Time-points named to reflect the average age of the children: 6mo, 6 months old; 18mo, 18 months old; 3yr, 3 years old; 5yr, 5 years old; 8yr, 8 years old; 14yr, 14 years old. Outcomes measured at all time-points are mother-rated, apart from at the 14yr time-point, where there are also child ratings, indicated by c.
Mean values:h2_vGWA = 0.02, Z = 1.04 h2_mGWA = 0.06, Z = 3.67 when Z>3:h2_vGWA = 0.06; h2_mGWA = * 0.08
datatable(
h2_TabFull %>% select(Trait,Age,Rater,I_mGWA,X2_mGWA,L_mGWA,I_vGWA,X2_vGWA,L_vGWA),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
Rel_h2 <-
ggplot(
h2_TabFull %>%
select(Trait, vGWA=h2_vGWA, mGWA=h2_mGWA, Reliability) %>%
pivot_longer(!c(Trait,Reliability),names_to="Effect",values_to="h2"),
aes(x=Reliability, y=h2, group=Effect, colour=Effect,shape=Effect,label=Trait)) +
geom_point() +
geom_smooth(method=lm,se=FALSE) +
scale_colour_viridis_d(option="viridis",direction=-1,begin=0.4,end=0.8)+
theme_minimal() +
directlabels::geom_dl(aes(label=Effect),method="smart.grid")+
guides(color="none",shape="none")#+
# ylab(bquote({h^2}))
plotly::ggplotly(Rel_h2, tooltip = c("Reliability","h2","Effect","Trait"))
S.Figure 3 shows the relationship between SNP-based heritability and scale reliability (Cronbach’s alpha) stratified by the type of analysis (GWA versus vGWA). Reliability, measured using Cronbach’s alpha, is the average correlation between each item and the total scale. For GWA, h2SNP increases with the reliability of the measure, whereas for vGWA, h2SNP decreases as the reliability of the measure increases.
mCorMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenFscAsr_rG_est.csv") %>% as.matrix()
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Ch","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q04","_6mo")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q05","_18mo")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q06","_3yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q07","_5yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Q09","_8yr")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat), "Q11C","_14yr(c)")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat), "Q11M","_14yr(m)")
row.names(mCorMat) <- colnames(mCorMat)
mCorTab <- round(mCorMat,2)
# mCorTab[upper.tri(mCorTab)] <- NA
mCorTab <- reshape2::melt(mCorTab)
names(mCorTab)[3] <- "rg"
ErrMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenFscAsr_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q04","_6mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q05","_18mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q06","_3yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q07","_5yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q09","_8yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11C","_14yr(c)")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11M","_14yr(m)")
row.names(ErrMat) <- colnames(ErrMat)
ErrTab <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"
mCorTab <- mCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))
mCorTab$Z <- mCorTab$rg/mCorTab$SE
Pval <- fread("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenFscAsr_rG_pval.txt")%>% select(Var1=V1,Var2=V2,Pval=V3) %>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q04","_6mo")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q05","_18mo")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q06","_3yr")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q07","_5yr")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q09","_8yr")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q11C","_14yr(c)")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q11M","_14yr(m)")) %>%
mutate(across(c(Var1,Var2),as.factor))
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
mCorTab <- mCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit() %>%
filter(!grepl(x=Var1,pattern=paste0(excl,collapse="|"))) %>%
filter(!grepl(x=Var2,pattern=paste0(excl,collapse="|"))) %>%
mutate(Z=round(Z,1))
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
mCor <- mCorTab %>% select(Var1,Var2,rg,SE,Z)
vars <- unique(c(levels(mCor$Var1),levels(mCor$Var2)))
excl <- c(
"IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
"CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
"SclFull","SppaFull","CbclAdhd"
)
vars <- vars[!grepl(paste0(excl,collapse="|"),vars)]
combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos <- data.frame(V1=vars,V2=vars)
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")
mCor <-
bind_rows(
mCorTab %>% select(Var1,Var2,rg),
mCorTab %>% select(Var1=Var2,Var2=Var1,rg)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(rg=ifelse(Var1==Var2,1,rg)) %>%
pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
column_to_rownames("Var2") %>%
as.matrix()
mCor[is.na(mCor)] <- 0
mCor[mCor>1] <- 1
mCor[mCor<(-1)] <- (-1)
Err <-
bind_rows(
mCorTab %>% select(Var1,Var2,SE),
mCorTab %>% select(Var1=Var2,Var2=Var1,SE)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(SE=ifelse(Var1==Var2,0,SE)) %>%
pivot_wider(names_from=Var1,values_from=SE) %>%
column_to_rownames("Var2") %>%
as.matrix()
Err[is.na(Err)] <- 0
Zsc <-
bind_rows(
mCorTab %>% select(Var1,Var2,Z),
mCorTab %>% select(Var1=Var2,Var2=Var1,Z)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(Z=ifelse(Var1==Var2,2,Z)) %>%
pivot_wider(names_from=Var1,values_from=Z) %>%
column_to_rownames("Var2") %>%
as.matrix()
Zsc[is.na(Zsc)] <- 0
Zsc <- abs(Zsc)
ChmCorplotly<-
heatmaply::heatmaply_cor(
mCor,
node_type = "scatter",
point_size_mat=Zsc,
point_size_name="|Z|",
label_names = c("x", "y", "rg"),
fontsize_col=8,
fontsize_row=8,
grid_size = 0.001,
dendrogram = "both",
subplot_heights=c(0.06, 0.94),
subplot_widths=c(0.94,0.06),
hide_colorbar=T
)
ChmCorplotly$height <- nrow(mCor)*14
ChmCorplotly$width <- nrow(mCor)*14
ChmCorplotly$x$data[[4]]$marker$size <- (ChmCorplotly$x$data[[4]]$marker$size/2.5)
ChmCorplotly
Genetic correlation matrices for mGWA (left) and vGWA (right) shows the degree to which genetic effects are shared across the phenotypes. The colour intensity represents the direction and strength of the genetic correlation (red = positive correlation, white = no correlation, blue = negative correlation) and the size of each point represents the z-statistic (rg / se). For mGWA (left), the correlation structure is more dense than for vGWA (right). This is likely because of differences in statistical power for detecting mean and variance effects - as a much smaller proportion of vGWA summary statistics had statistically significant h2 estimates.
vCorMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenQrsAsr_rG_est.csv") %>% as.matrix()
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Ch","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q04","6mo")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q05","18mo")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q06","3yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q07","5yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q09","8yr")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q11C","14yr(c)")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Q11M","14yr(m)")
row.names(vCorMat) <- colnames(vCorMat)
vCorTab <- round(vCorMat,2)
# vCorTab[upper.tri(vCorTab)] <- NA
vCorTab <- reshape2::melt(vCorTab)
names(vCorTab)[3] <- "rg"
ErrMat <- fread("Tables_GenSem/ChPsy/MvLDSC_ChPsyGenQrsAsr_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q04","6mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q05","18mo")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q06","3yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q07","5yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q09","8yr")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11C","14yr(c)")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Q11M","14yr(m)")
row.names(ErrMat) <- colnames(ErrMat)
ErrTab <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"
vCorTab <- vCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))
vCorTab$Z <- vCorTab$rg/vCorTab$SE
Pval <- fread("Tables_GenSem/ChPsy/mvLDSC_ChPsyGenQrsAsr_rG_pval.txt")%>%
select(Var1=V1,Var2=V2,Pval=V3) %>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),as.factor)) %>%
mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q04","6mo")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q05","18mo")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q06","3yr")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q07","5yr")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q09","8yr")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q11C","14yr(c)")) %>%
mutate(across(c(Var1,Var2),str_replace_all,"Q11M","14yr(m)")) %>%
mutate(across(c(Var1,Var2),as.factor))
vCorTab <- vCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()%>%
mutate(across(c(Var1,Var2),as.character))%>%
filter(!grepl(x=Var1,pattern=paste0(excl,collapse="|"))) %>%
filter(!grepl(x=Var2,pattern=paste0(excl,collapse="|"))) %>%
mutate(Z=round(Z,1))
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vCor <- vCorTab %>% select(Var1,Var2,rg,SE,Z)
vars <- unique(c(vCor$Var1,vCor$Var2))
excl <- c(
"IcqFull","CbclFull","EasFull","CfqFull","ItseaFull","ScqFull","SdqFull",
"CccFull","SprakFull","NhpicFull","RsdbdFull","PeqFull","SbFull",
"SclFull","SppaFull","CbclAdhd"
)
vars <- vars[!grepl(paste0(excl,collapse="|"),vars)]
combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos <- data.frame(V1=vars,V2=vars)
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")
vCor <-
bind_rows(
vCorTab %>% select(Var1,Var2,rg),
vCorTab %>% select(Var1=Var2,Var2=Var1,rg)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(rg=ifelse(Var1==Var2,1,rg)) %>%
pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
column_to_rownames("Var2") %>%
as.matrix()
vCor[is.na(vCor)] <- 0
vCor[vCor>1] <- 1
vCor[vCor<(-1)] <- (-1)
Err <-
bind_rows(
vCorTab %>% select(Var1,Var2,SE),
vCorTab %>% select(Var1=Var2,Var2=Var1,SE)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(SE=ifelse(Var1==Var2,0,SE)) %>%
pivot_wider(names_from=Var1,values_from=SE) %>%
column_to_rownames("Var2") %>%
as.matrix()
Err[is.na(Err)] <- 0
Zsc <-
bind_rows(
vCorTab %>% select(Var1,Var2,Z),
vCorTab %>% select(Var1=Var2,Var2=Var1,Z)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(Z=ifelse(Var1==Var2,2,Z)) %>%
pivot_wider(names_from=Var1,values_from=Z) %>%
column_to_rownames("Var2") %>%
as.matrix()
Zsc[is.na(Zsc)] <- 0
Zsc <- abs(Zsc)
ChvCorplotly<-
heatmaply::heatmaply_cor(
as.data.frame(vCor),
node_type = "scatter",
# point_size_mat = -log(Pval),
# point_size_name = "-logP",
point_size_mat=Zsc,
point_size_name="|Z|",
label_names = c("x", "y", "rg"),
fontsize_col=8,
fontsize_row=8,
grid_size = 0.001,
dendrogram = "both",
subplot_heights=c(0.06, 0.94),
subplot_widths=c(0.94,0.06),
hide_colorbar=T
)
ChvCorplotly$height <- nrow(vCor)*17
ChvCorplotly$width <- nrow(vCor)*17
ChvCorplotly$x$data[[4]]$marker$size <- (ChvCorplotly$x$data[[4]]$marker$size/2)
ChvCorplotly
Genetic correlation matrices for mGWA (left) and vGWA (right) shows the degree to which genetic effects are shared across the phenotypes. The colour intensity represents the direction and strength of the genetic correlation (red = positive correlation, white = no correlation, blue = negative correlation) and the size of each point represents the z-statistic (rg / se). For mGWA (left), the correlation structure is more dense than for vGWA (right). This is likely because of differences in statistical power for detecting mean and variance effects - as a much smaller proportion of vGWA summary statistics had statistically significant h2 estimates.
RgM_N <- nrow(mCorTab %>% filter(Z > 3))
meanMRgZ <- round(mean(abs(mCorTab$Z)),2)
datatable(
mCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
filter='bottom',rownames=F,extensions='Buttons',
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
meanVRgZ <- round(mean(abs(vCorTab$Z)),2)
RgV_N <- nrow(vCorTab %>% filter(Z > 3))
datatable(
vCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
SnpTab <- fread("Tables_Snps/ChAnth/TopSnps_Clumped.tsv") %>%
mutate(af = round(af,2)) %>%
mutate(mP = formatC(mP, format = "e", digits = 2)) %>%
mutate(vP = formatC(vP, format = "e", digits = 2)) %>%
mutate(across(c(Trait),str_remove_all,"Ch"))
DT::datatable(
SnpTab,
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
CorP <- cor.test(as.numeric(SnpTab$mP),as.numeric(SnpTab$vP))
Correlation between p-values for mGWA and vGWA: Pearson’s r = 2.68; p = 7.50e-03
vPheWasTab <- fread("Tables_Snps/ChAnth/Anth_vPheWasTab.tsv")
mPheWasTab <- fread("Tables_Snps/ChAnth/Anth_mPheWasTab.tsv")
PheWasTab3 <- vPheWasTab %>% bind_rows(mPheWasTab)
PheWasLab <-
PheWasTab3 %>%
arrange(PheWasP) %>% distinct(PheWasTrait,MoBaTrait,Age,.keep_all=T) %>%
select(PheWasTrait,MoBaTrait,Age,PheWasP) %>%
mutate(Label=str_replace_all(PheWasTrait, "\\s*\\([^\\)]+\\)", "")) %>%
mutate(Label=str_replace_all(Label, "[[:punct:]]", "")) %>%
mutate(Label=str_to_title(Label)) %>%
mutate(Label=str_replace_all(Label, " ", "")) %>%
mutate(Label=str_replace_all(Label, "CellCount", "Count")) %>%
# distinct(Label,MoBaTrait,.keep_all=T) %>%
inner_join(PheWasTab3,c("MoBaTrait","PheWasTrait","PheWasP","Age")) %>%
select(chr,rsid,MoBaTrait,Age,PheWasTrait=Label,PheWasP,vP,mP) %>% distinct() %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "C_", " ")) %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "Full", "")) %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([A-Z])", "\\1 \\2")) %>%
mutate(MoBaTrait=str_replace_all(MoBaTrait, "([a-z])([0-9])", "\\1 \\2"))%>%
mutate(Age = factor(Age, levels=c("18mo","03yo","05yo","07yo","08yo","14yo")))
PheWasLab <-
bind_rows(
PheWasLab %>% filter(vP<1e-5) %>% rename("GwaP"=vP) %>% mutate(Effect="vQTL"),
PheWasLab %>% filter(mP<1e-5) %>% rename("GwaP"=mP) %>% mutate(Effect="mQTL")
) %>%
mutate(GwaP=round(-log10(GwaP),1),PheWasP=round(-log10(PheWasP),1))
PheWasLab$mP <- formatC(PheWasLab$mP, format = "e", digits = 2)
PheWasLab$vP <- formatC(PheWasLab$vP, format = "e", digits = 2)
PheWasPlot <-
ggplot(PheWasLab,
aes(y=PheWasP,x=GwaP,shape=Effect,colour=MoBaTrait,label=PheWasTrait))+
geom_point()+theme_minimal()+
facet_grid(cols=vars(Age),rows=vars(Effect),drop=T)+
scale_colour_viridis_d(begin=.2,end=.8)+
scale_shape_manual(values=c(0:14))+
theme(legend.position="none")
ggplotly(PheWasPlot)
DT::datatable(
PheWasLab %>% select(chr,rsid,MoBaTrait,Age,PheWasTrait,PheWasP,vP,mP),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX=T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
h2TabM <-
fread(paste0("Tables_GenSem/ChAnth/mvLDSC_ChAnth_h2FullTable.txt"),header=T) %>%
mutate(across(c(h2, SE, Z), ~round(.x,2))) %>%
distinct(Trait,.keep_all = T) %>%
separate(Trait,c("Trait","Age","Adj"),"_") %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
filter(!(Trait=="Weight"&is.na(Adj))) %>%
mutate(Type = paste("mGWA")) %>%
select(-Adj)
h2TabV <-
fread(paste0("Tables_GenSem/ChAnth/mvLDSC_ChAnthQrs_h2FullTable.txt"),header=T) %>%
mutate(across(c(h2, SE, Z), ~round(.x,2))) %>%
distinct(Trait,.keep_all = T) %>%
separate(Trait,c("Trait","Age","Adj"),"_") %>%
mutate(Trait=str_replace_all(Trait,"Ch","")) %>%
filter(!(Trait=="Weight"&is.na(Adj))) %>%
mutate(Type = paste("vGWA")) %>%
select(-Adj)
h2Tab1 <- rbind(h2TabM,h2TabV)
h2Tab2 <-
h2TabM %>%
select(Trait,Age,h2_mGWA=h2,SE_mGWA=SE,Z_mGWA=Z) %>%
full_join(h2TabV %>% select(Trait,Age,h2_vGWA=h2,SE_vGWA=SE,Z_vGWA=Z))
N <- fread("./Tables_GenSem/ChAnth/ChAnthQrs_SampleSizes.txt",header=F) %>%
mutate(V1=str_replace_all(V1,"_RegenieMQt.Munged.sumstats.gz","")) %>%
mutate(V1=str_replace_all(V1,"_RegenieVQt.Munged.sumstats.gz","")) %>%
mutate(V1=str_replace_all(V1,"GwaSumStats/ChAnthQrs/","")) %>%
# mutate(V1=str_replace_all(V1,"([a-z])([A-Z])", "\\1 \\2")) %>%
separate(V1, c("Trait","Age","Adj")) %>%
select(Trait, Age, N=V2) %>%
mutate(Trait=str_replace_all(Trait,"Ch",""))
meanh2V <- round(mean(h2TabV$h2),2)
meanh2m <- round(mean(h2TabM$h2),2)
meanZV <- round(mean(h2TabV$Z),2)
meanZm <- round(mean(h2TabM$Z),2)
meanh2Vz <- round(mean(h2TabV %>% filter(Z>3) %>% .$h2),2)
meanh2mz <- round(mean(h2TabM %>% filter(Z>3) %>% .$h2),2)
h2Tab2$Age <- factor(h2Tab2$Age,levels=c("18mo","03yo","05yo","07yo","08yo","14yo"))
datatable(
h2Tab2, # %>% mutate(`|Zm|` = abs()),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
meanh2V <- round(mean(h2TabV$h2),4)
maxh2V <- round(max(h2TabV$h2),4)
meanh2m <- round(mean(h2TabM$h2),4)
maxh2m <- round(max(h2TabM$h2),4)
meanZV <- round(mean(h2TabV$Z),2)
meanZm <- round(mean(h2TabM$Z),2)
meanh2Vz <- round(mean(h2TabV %>% filter(Z>3) %>% .$h2),2)
meanh2mz <- round(mean(h2TabM %>% filter(Z>3) %>% .$h2),2)
h2V_N <- nrow(h2TabV %>% filter(Z > 3))
h2M_N <- nrow(h2TabM %>% filter(Z > 3))
meanXV <- round(mean(h2TabV$MeanChi2),2)
meanXm <- round(mean(h2TabM$MeanChi2),2)
Mean h2-vGWA = 0.0208, meanZ= 1.09; Mean
h2-mGWA = 0.3008, meanZ= 10.11;
when Z>3: Mean h2-vGWA = NaN; Mean h2-mGWA =
0.32
Mh2_vGWAh2 <-
ggplot(
h2Tab2 %>%
mutate(
h2_GWA=round(h2_mGWA,2),h2_vGWA=round(h2_vGWA,2),
se_GWA=round(SE_mGWA,2),se_vGWA=round(SE_vGWA,2),
lse_vGWA=round(h2_vGWA-SE_vGWA,2),use_vGWA=round(h2_vGWA+SE_vGWA,2),
lse_GWA=round(h2_mGWA-SE_mGWA,2), use_GWA=round(h2_mGWA+SE_mGWA,2)
),
aes(
x=h2_GWA,y=h2_vGWA
)) +
geom_linerange(aes(ymin=round(h2_vGWA-SE_vGWA,2),ymax=round(h2_vGWA+SE_vGWA,2),color=Trait)) +
geom_linerange(aes(xmin=round(h2_mGWA-SE_mGWA,2),xmax=round(h2_mGWA+SE_mGWA,2),color=Trait)) +
scale_x_continuous(limits=c(-0.1,0.5),breaks=c(seq(0,0.5,0.1)))+
scale_y_continuous(limits=c(-0.1,0.5),breaks=c(seq(0,0.5,0.1)))+
geom_hline(yintercept=0,colour="limegreen",linetype="dotted") +
geom_vline(xintercept=0,colour="limegreen",linetype="dotted") +
# geom_smooth(method=lm,se=FALSE,colour="limegreen") +
facet_wrap(facets=vars(Age),ncol=3) + #, scales = "free",space = "free")+
theme_minimal() +
scale_colour_viridis_d(option="plasma",begin=0.2,end=0.8)+
guides(color="none")
ggplotly(Mh2_vGWAh2, tooltip = c("Trait","h2_GWA","se_GWA","h2_vGWA","se_vGWA"))
mCorMat <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnth_rG_est.csv") %>% as.matrix()
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(mCorMat) <- str_replace_all(colnames(mCorMat),"Ch","")
row.names(mCorMat) <- colnames(mCorMat)
KeepRows <- which(grepl(pattern="Height|Adj",rownames(mCorMat)))
KeepCols <- which(grepl(pattern="Height|Adj",colnames(mCorMat)))
mCorMat <- mCorMat[KeepRows,KeepCols]
colnames(mCorMat) <- str_remove_all(string=colnames(mCorMat),pattern="_hAdj")
row.names(mCorMat) <- colnames(mCorMat)
mCorTab <- round(mCorMat,2)
# mCorTab[upper.tri(mCorTab)] <- NA
mCorTab <- reshape2::melt(mCorTab)
names(mCorTab)[3] <- "rg"
ErrMat <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnth_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
row.names(ErrMat) <- colnames(ErrMat)
ErrMat <- ErrMat[KeepRows,KeepCols]
colnames(ErrMat) <- str_remove_all(string=colnames(ErrMat),pattern="_hAdj")
row.names(ErrMat) <- colnames(ErrMat)
ErrTab <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"
mCorTab <- mCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))
mCorTab$Z <- mCorTab$rg/mCorTab$SE
Pval <- fread("Tables_GenSem/ChAnth/mvLDSC_ChAnth_rG_pval.txt")%>% select(Var1=V1,Var2=V2,Pval=V3) %>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>%
mutate(across(c(Var1,Var2),as.factor)) %>%
filter(grepl(pattern="Height|Adj",Var1))%>%
filter(grepl(pattern="Height|Adj",Var2)) %>%
mutate(across(c(Var1,Var2),str_remove_all,"_hAdj"))
mCorTab <- mCorTab %>%
full_join(Pval,by=c("Var1","Var2")) %>% na.omit() %>%
mutate(Z=round(Z,1))
mCor <- mCorTab %>% select(Var1,Var2,rg,SE,Z)
vars <- unique(c(mCor$Var1,mCor$Var2))
combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos <- data.frame(V1=vars,V2=vars)
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")
mCor <-
bind_rows(
mCorTab %>% select(Var1,Var2,rg),
mCorTab %>% select(Var1=Var2,Var2=Var1,rg)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(rg=ifelse(Var1==Var2,1,rg)) %>%
pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
column_to_rownames("Var2") %>%
as.matrix()
mCor[is.na(mCor)] <- 0
mCor[mCor>1] <- 1
mCor[mCor<(-1)] <- (-1)
Err <-
bind_rows(
mCorTab %>% select(Var1,Var2,SE),
mCorTab %>% select(Var1=Var2,Var2=Var1,SE)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(SE=ifelse(Var1==Var2,0,SE)) %>%
pivot_wider(names_from=Var1,values_from=SE) %>%
column_to_rownames("Var2") %>%
as.matrix()
Err[is.na(Err)] <- 0
Zsc <-
bind_rows(
mCorTab %>% select(Var1,Var2,Z),
mCorTab %>% select(Var1=Var2,Var2=Var1,Z)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(Z=ifelse(Var1==Var2,2,Z)) %>%
pivot_wider(names_from=Var1,values_from=Z) %>%
column_to_rownames("Var2") %>%
as.matrix()
Zsc[is.na(Zsc)] <- 0
Zsc <- abs(Zsc)
ChmCorplotly<-
heatmaply::heatmaply_cor(
mCor,
node_type = "scatter",
# point_size_mat = -log(Pval),
# point_size_name = "-logP",
point_size_mat=Zsc,
point_size_name="|Z|",
label_names = c("x", "y", "rg"),
fontsize_col=8,
fontsize_row=8,
grid_size = 0.001,
dendrogram = "both",
subplot_heights=c(0.06, 0.94),
subplot_widths=c(0.94,0.06),
hide_colorbar=T
)
# ChmCorplotly$height <- nrow(mCor)*15
# ChmCorplotly$width <- nrow(mCor)*15
# ChmCorplotly$x$data[[4]]$marker$size <- (ChmCorplotly$x$data[[4]]$marker$size/2)
ChmCorplotly
vCorMat <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnthQrs_rG_est.csv") %>% as.matrix()
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(vCorMat) <- str_replace_all(colnames(vCorMat),"Ch","")
row.names(vCorMat) <- colnames(vCorMat)
KeepRows <- which(grepl(pattern="Height|Adj",rownames(vCorMat)))
KeepCols <- which(grepl(pattern="Height|Adj",colnames(vCorMat)))
vCorMat <- vCorMat[KeepRows,KeepCols]
colnames(vCorMat) <- str_remove_all(string=colnames(vCorMat),pattern="_hAdj")
row.names(vCorMat) <- colnames(vCorMat)
vCorTab <- round(vCorMat,2)
# vCorTab[upper.tri(vCorTab)] <- NA
vCorTab <- reshape2::melt(vCorTab)
names(vCorTab)[3] <- "rg"
ErrMat <- fread("Tables_GenSem/ChAnth/MvLDSC_ChAnthQrs_rG_err.csv") %>% as.matrix()
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieMQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"_RegenieVQt.Munged.sumstats.gz","")
colnames(ErrMat) <- str_replace_all(colnames(ErrMat),"Ch","")
row.names(ErrMat) <- colnames(ErrMat)
ErrMat <- ErrMat[KeepRows,KeepCols]
colnames(ErrMat) <- str_remove_all(string=colnames(ErrMat),pattern="_hAdj")
row.names(ErrMat) <- colnames(ErrMat)
ErrTab <- round(ErrMat,3)
# ErrTab[upper.tri(ErrTab)] <- NA
ErrTab <- reshape2::melt(ErrTab)
names(ErrTab)[3] <- "SE"
vCorTab <- vCorTab %>% full_join(ErrTab,by=c("Var1","Var2"))
vCorTab$Z <- vCorTab$rg/vCorTab$SE
Pval <- fread("Tables_GenSem/ChAnth/mvLDSC_ChAnthQrs_rG_pval.txt")%>%
select(Var1=V1,Var2=V2,Pval=V3) %>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieVQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),str_replace_all,"_RegenieMQt.Munged.sumstats.gz",""))%>%
mutate(across(c(Var1,Var2),as.factor)) %>%
mutate(across(c(Var1,Var2),str_remove_all,"Ch")) %>%
mutate(across(c(Var1,Var2),as.factor)) %>%
filter(grepl(pattern="Height|Adj",Var1))%>%
filter(grepl(pattern="Height|Adj",Var2)) %>%
mutate(across(c(Var1,Var2),str_remove_all,"_hAdj"))
vCorTab <- vCorTab %>% full_join(Pval,by=c("Var1","Var2")) %>% na.omit()%>%
mutate(across(c(Var1,Var2),as.character))%>%
mutate(Z=round(Z,1))
vCor <- vCorTab %>% select(Var1,Var2,rg,SE,Z)
vars <- unique(c(vCor$Var1,vCor$Var2))
combos <- apply(combn(vars,2),1,paste) %>% as.data.frame()
autos <- data.frame(V1=vars,V2=vars)
combos <- combos %>% bind_rows(autos) %>% arrange(V1,V2)
names(combos)<-c("Var1","Var2")
vCor <-
bind_rows(
vCorTab %>% select(Var1,Var2,rg),
vCorTab %>% select(Var1=Var2,Var2=Var1,rg)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(rg=ifelse(Var1==Var2,1,rg)) %>%
pivot_wider(names_from=Var1,values_from=rg,values_fill=NA) %>%
column_to_rownames("Var2") %>%
as.matrix()
vCor[is.na(vCor)] <- 0
vCor[vCor>1] <- 1
vCor[vCor<(-1)] <- (-1)
Err <-
bind_rows(
vCorTab %>% select(Var1,Var2,SE),
vCorTab %>% select(Var1=Var2,Var2=Var1,SE)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(SE=ifelse(Var1==Var2,0,SE)) %>%
pivot_wider(names_from=Var1,values_from=SE) %>%
column_to_rownames("Var2") %>%
as.matrix()
Err[is.na(Err)] <- 0
Zsc <-
bind_rows(
vCorTab %>% select(Var1,Var2,Z),
vCorTab %>% select(Var1=Var2,Var2=Var1,Z)
) %>%
full_join(combos,c("Var1","Var2"))%>% arrange(Var1,Var2) %>%
mutate(Z=ifelse(Var1==Var2,2,Z)) %>%
pivot_wider(names_from=Var1,values_from=Z) %>%
column_to_rownames("Var2") %>%
as.matrix()
Zsc[is.na(Zsc)] <- 0
Zsc <- abs(Zsc)
ChvCorplotly<-
heatmaply::heatmaply_cor(
as.data.frame(vCor),
node_type = "scatter",
# point_size_mat = -log(Pval),
# point_size_name = "-logP",
point_size_mat=Zsc,
point_size_name="|Z|",
label_names = c("x", "y", "rg"),
fontsize_col=8,
fontsize_row=8,
grid_size = 0.001,
dendrogram = "both",
subplot_heights=c(0.06, 0.94),
subplot_widths=c(0.94,0.06),
hide_colorbar=T
)
ChvCorplotly
datatable(
mCorTab %>% select(Var1,Var2,rg,SE,Pval,Z) ,
filter='bottom',rownames=F,extensions='Buttons',
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))
datatable(
vCorTab %>% select(Var1,Var2,rg,SE,Pval,Z),
filter='bottom',extensions='Buttons',rownames=F,
options=list(
dom='lfrtiBp',buttons=c('csv'),
searching=F,pageLength=10,lengthMenu=list(c(10,-1),c('10','All')),scrollX = T,
columnDefs=list(list(className='dt-left',targets="_all"))
))